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Abstract. In this work, the quantum phase transition in the sub-Ohmic spin-boson model is studied using a single- 
mode approximation, by combining the rotating wave transformation and the transformations used in the numerical 
renormalization group (NRG). Analytical results for the critical coupling strength as c , the magnetic susceptibility 
and the spin-spin correlation function C(w) at finite temperatures are obtained and further confirmed by numerical 
results. We obtain the same a c as the mean-field approximation. The critical exponents are classical: f} = 1/2, 6 = 3, 
7 = 1, x = 1/2, y* = 1/2, in agreement with the spin-boson model in < s < 1/2 regime. C(u) has nontrivial behavior 
reflecting coherent oscillation with temperature dependent damping effects due to the environment. We point out the 
original NRG has problem with the crossover temperature T*, and propose a chain Hamiltonian possibly suitable for 
implementing NRG without boson state truncation error. 

PACS. 05.30.Jp-05.10.Cc-64.70.Tg 



1 Introduction 

The spin-boson model (SBM) is the simplest model that de- 
scribes a quantum two-level system coupled to a dissipative en- 
vironment. SBM has extensive relevance to physical systems in 
condensed matter physics QUI, quantum optics, and quantum 
chemistry. Its properties are studied extensively. For the bath 
spectra exponent < s < 1 (sub-ohmic bath), the ground state 
may change from the spin-tunneling state to the spin-pinned 
state through a second order phase transition, as the dissipa- 
tion strength increase above a critical value. This environment- 
induced quantum phase transition attracts much attention in the 
past few years 



and the universality class 
of this transition is under intensive studies I113II14|[T51|16II17I 
[TM9ll20ll2T1l22ll23ll. 

The standard quantum-classical mapping predicts that the 
sub-Ohmic SBM is equivalent to a classical Ising spin chain 
with long-range interaction [1]. Classical Monte Carlo simu- 
lations 11241 and renormalization group analysis Il25ll for the 
long-range Ising model predict a continuous magnetic transi- 
tion with classical critical exponents for < s < 1/2, i.e., 
P — 1/2, 6 = 3, and y = 1. In Ref. [4j, using the numerical 
renormalization group (NRG) technique, a continuous transi- 
tion was found in the sub-Ohmic SBM, with non-classical ex- 
ponents for < s < 1/2. Recently, other methods such as 
quantum monte carlo (QMC) [16], sparse polynomial space ap- 
proach (8), and extended coherent state techniques [6,7] have 
found classical critical exponents for < s < 1/2. 

Later research disclosed that although non-universal quan- 
tities such as a c are obtained quite accurately, NRG results are 
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not reliable for critical exponents ft, 6, and x in the deep sub- 
Ohmic regime < s < 1/2. This is due to two inherent limi- 
tations of the bosonic NRG (BNRG): the boson Hilbert-space 
truncation error and the mass flow error fl5l . On one hand, 
the BNRG employs boson Hilbert-space truncation in the iter- 
ative diagonalization process. It causes errors in ft and 6 which 
characterize the flow into the localized phase at zero temper- 
ature M15II23L On the other hand, the mass flow causes errors 
for non-zero temperatures as a result of neglecting of low-lying 
bath modes with energy smaller than temperature 11811191 . 

Recently, in Ref. l23l . the BNRG calculation is combined 
with finite size scaling analysis of Nb, the number of truncated 
boson states, to recover the correct exponents ft and 6. Using 
the matrix product state to describe the ground state of the Wil- 
son chain, Guo et al. were able to take the optimal boson ba- 
sis and successfully overcome the boson state truncation error. 
However, for accurate study of finite temperature as well as dy- 
namical properties, a NRG-like algorithm without boson state 
truncation is still desirable, with its energy flow incorporating 
physical information of excited states. As an attempt to find 
possible NRG variants for this purpose, we carry out the rotat- 
ing wave transformation (RWT) II 1011 11 111 211261 to the Hamilto- 
nian of the SBM, to single out the boson displacement induced 
by the spin-boson coupling, so that the truncation errors for the 
displaced boson modes could be avoided in the NRG calcu- 
lation. In the limit of complete transformation r] — > (rj is a 
control parameter of RWT), a boson mode is singled out in the 
Wilson chain Hamiltonian, which has vanishing frequency and 
hopping amplitude to other sites. In this paper, we consider the 
single mode Hamiltonian H^irf) in which the spin is coupled 
only to this adiabatic mode, as the first step towards the BNRG 
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study of the full chain. Using adiabatic approximation which 
is exact in the adiabatic limit 77 — > and checking it using ex- 
act diagonalization method, we find that this model gives the 
critical coupling strength a c - sA/2u> c , which is same as the 
mean-field result 11231 . The critical exponents are obtained as 
/3 — 1/2, 6 = 3, y = 1, x = 1/2, and y* = 1/2, in agreement 
with the expected classical exponents of the SBM in the regime 
< s < 1/2. 

The remainder of the paper is organized as follows: In Sec. 2, 
we introduce the SBM and derive the approximate single-mode 
Hamiltonian Ho(r]) for it. In Section 3, we present results for the 
order parameter (<r z ), susceptibility %(T) and spin-spin corre- 
lation function C(lo). Related issues including the scaling re- 
lations, the BNRG result T*, and a possible NRG algorithm 
without boson state truncation error are discussed. Finally, we 
close with a summary. 



where x = IZkgk {a\ - a k j. 

Similar RWT's have been widely used in the study of SBM. 
Combined with perturbation 01 1 11261 or numerical diagonal- 
ization J6), both critical coupling strength a c and dynamical 
properties are accurately obtained, 77 measures how much dis- 
placement in the k-th mode induced by coupling to cr z is ab- 
sorbed into the new boson mode. For 77 = 0, the coupling term 
o~ z (a\ + ak) disappears in H. This term is the source of the dis- 
placement divergence in low energy boson modes, when the 
coupling is strong. It is expected that in the limit 77 — > 0, RWT 
removes most of the boson displacement and it thus helps to 
avoid the boson state truncation error in BNRG. Indeed, for the 
case A — 0, H(t] = 0) has no displaced bosons. 

For H, we define the bath spectral function as 



2 Model and Method 

2.1 effective single mode Hamiltonian H Q (r/) 
The Hamiltonian of SBM reads(/? = 1) HE 

A e v— 1 j. <x, v— ' t 
H S b = - 2 cr x + -cr z + ^ cotaiai + — ^ Ma- + ad. (1) 



Here, the two-state system is represented by the Pauli matrices 
cr x and <t z . e is the bias of the two states and A is the tunneling 
strength between them. The environment is modeled by a col- 
lection of harmonic oscillators, with creation and annihilation 
operators a. and a,, respectively, w, and /I, represent the fre- 
quency of the i-th boson mode and its coupling to the two-state 
system, respectively. The effect of the harmonic environment is 
characterized by the bath spectral function 



J(oj) = n ^ Aj6(u> - a>i). 



(2) 



We use a model spectral function as 



J(aS) = 2nau) l c "cj s , (0 < u < a) c , s > -1), (3) 

where the dimensionless parameter a characterizes the dissipa- 
tion strength, and the cutoff <d c = 1 is set as the energy unit. 

Starting from Eq.lQ), we carry out RWT fl0l[TTi nll26l H = 
UHsbU' 1 . Here U = e s and § = l/2£jt£fc(<z! - a k )cr z . g k is 
the parameter controlling the transformation for the fe-th boson 
mode, for which we use 



gk = 



OJ k +T] 



(4) 



Here 77 is the parameter to control the RWT. ft is then written 
as (neglecting a constant) 

H = -cr z + 2^ u k a\a k + -cr- ^ gk(a[ + a k ) 



A i 
- -cr x cosh (x) - -Ao-y sinh (x), 



(5) 



J(u). 



(6) 



Note that in Eq.©, the local boson modes that are coupled to 
<t z and to <x v are the same, i.e., £t g k a k - This is a result of the 
special form of g k in Eq.©. With possible NRG-related appli- 
cations in mind, we further carry out logarithmic discretization 
for Eq.©, following the standard procedure [4 j. We obtain the 
star-Hamiltonian as 



H, = -a, + 



+00 ^ +00 

V £„ala„ + — —cr z V y n (a„ + a\ 



A iA 

-a x cosh(£) -y(Tj sinhCf-), 



(7) 



with 



-If ' 

n — j I 

r>A~ n ti) c 



J{oS)iodcj 



J(u)du>. 



(8) 



In Eq.dTJi, ^ = 1/(77 V^) IiZ=o ln{a l n - a„). A > 1 is the logarith- 
mic parameter. 

In Fig.l, we plot J((l>) on log scale and the coefficients g n 
and y n jr\ for various rfs at s = 0.3. J(a>) is characterized by a 
crossover scale of ~ 77, which separates J(ai) ~ a> s in a> « to* 



and J(oj) 



in to » of regime. Correspondingly, y n 



y^(«/2)(i-s) m sma u n anc [ y n ^ ^-(n/2)(i+*) m laj-gg n regime. is 
independent of 77 in the limit 77 — > 0. In the limit 77 — > 00, J(to) = 
J((l>) and H s recovers the original star-Hamiltonian of SBM [4 1 . 
In the limit 77 — > 0, J(oS) — > 0. y n disappears but y„/rj follows 
y^(«/2)(i-s) an( j di ver g es j n large 77 limit. As a result, the coupling 

term y„cr z (a„ + a^) in H s will disappear as 77 — > 0. The role of 
coupling operators are now played by cosh(£) and sinh(£) in 
which x is a local boson mode with diverging coefficient. 

Using orthogonal transformation, the star-Hamiltonian can 
be transformed into a semi-infinite chain which is suitable for 
iterative diagonalization in NRG ||29l . 

He = I cr~ + - J— o", (b + b' ] ) - i-a- x cosh (x) 
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Fig. 1. (a) Transformed boson spectral function J(uj). From top to 
bottom: 77 = 10 2 , 10 1 , 10°, 10 _1 , 10~ 2 , 1CT 3 , 10~ 4 . (b) £„ (dashed line 
independent of 77) and y„/r] (solid lines). From bottom to top: 77 = 
10 2 , 10°, lfr 2 , lfr\ lO" 6 . Other parameters are 5 = 03, A = OA, a = 
0.03, and A = 2.0 



- -icr y sinhCf) + ^ \s n b\b„ + t n {V n b n + x + ^ +1 &«)] .(9) 

n=0 

with;? = I h V^o/tt (b l - b ). The parameters 770 = Yjnjl- e« 
and t n 's are obtained through the recursion relations as |4| 



n=0 



n=0 

Urn+ln ~ — £m)Umn ~ hn-\^ m-\n\ • 



(10) 



The initial conditions are f_i = 0, I/_i„ = 0, and Uq„ — y n j y[rjo~. 

Now we analyze the asymptotic form of the coefficients in 
H c . The coupling constant behaves as 770 oc tj 1+s so that the term 
cr z (b Q + /bo) disappears in the limit 77 — > 0. The strong coupling 
fixed point at A = can then be obtained correctly because 
no displacements of bosons are present. The price is that the 
parameter in;? diverges as sfrjo/ri oc rf~ l+s V 2 . In Fig. 2(a), we 
plot e n A n and t„A" as functions of n, for various 77's at s — 0.3. 
It is seen that e„ and t„ have a plateau in high energy regime 
n < n cr . While for n > n cr , they approach the levels same as in 
the original NRG chain. Here, n cr is a crossover scale related 
to the of of J(a>). The plateaus in the curves extend to larger 
n for smaller 77. In the limit 77 — > 0, they form new levels. This 
means that e„(77 = 0) and t n (r\ = 0) have different values from 
e n {r\ = 00) and t n (r\ - 00) but still decays as A~". This is a 
necessary condition for NRG to be applicable for H c (j] = 0). 

Direct NRG calculation based on H c shows that although 
the original BNRG results is recovered using 77 > 100, the flow 
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Fig. 2. (a) Chain coefficients e„A" (solid lines) and t„A" (dashed lines) 
as functions of 77 for various 77's. From bottom to top, 77 = 10 1 , 10°, 
1CT 7 . (b) Some coefficients as functions of 77: 770 (squares), eo (cir- 
cles), 61 (up triangles), to (down triangles), and fi (diamonds). Other 
parameters are S = 0.3, A = 0.1, a = 0.03, and A = 2.0 



lines fails to describe a fixed point for 77 < 10 2 . The reason 
lies in the following asymptotic relations as shown in Fig. 2(b), 



eo oc 77 





t() oc 77 

,0 



(l-j)/2 



e\ oc 77 , t\ oc 77 



(11) 



In the limit 77 = 0, eo = to = 0. This means that the zero- 
th mode of the boson chain becomes an adiabatic boson mode 
with A/eo — > 00. At the same time, it tends to get disconnected 
to the rest of the chain due to fo — > 0. The Hamiltonian Ho 
which contains the spin and the 0-th boson mode reads 



#0 = ^cr ; + \ J—cr z (bo + b Q ) + sob J bo 



A A 
- -<r x cosh (£) - -10-y sinh (£). 



(12) 



When doing iterative diagonalization within NRG, Ho is first 
diagonalized and the second boson mode b\ is added to form 
a new Hamiltonian. In the limit e — 0, the energy levels of 
Ho tends to be infinitely dense and is close to degeneracy. As a 
result, the next boson site to be added will not be a small pertur- 
bation to Hq. The basis for NRG to work, i.e., the exponential 
decay of energy scale, is not fulfilled. 

Although the most direct endeavor of taking 77 — » does 
not give a useful chain Hamiltonian for BNRG, the obtained 
single-mode Hamiltonian Ho is still interesting. First, it repre- 
sents the limit of long-range interaction in the imaginary time 
axis, and its critical properties should show the characteristics 
of the SBM in small s limit. Second, Hq describes a spin cou- 
pled to an adiabatic resonator which by itself is an interesting 
system. Previous study in terms of the entanglement entropy 
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shows that a quantum phase transition occurs in this system 
|27|. Third, as an exactly solvable model, Ho serves as the zero- 
th order approximation to H c in the sense that the full model 
can be studied by treating the hopping term fo {b Q b\ + b\bo) as 
a small perturbation. Therefore, in the following sections, we 
will focus on the properties of Ho(t] — > 0) and leave the full 
Hamiltonian H c for a future study. 



For numerical calculations, we start from Ho in Eq. (fT2l) . We 
build matrices for cosh(£ ) and sinh(£ ) on Nb boson occupation 
bases {|0 >, |1 >, \Nb - 1 >}, using the eigen wave func- 
tions of a displaced harmonic oscillator 11281 . In the adiabatic 
limit, the displacement is strong. We use an efficient algorithm 
to produce the associated Laguerre polynomial for calculating 
the wave functions. The details are presented in Appendix B. 



2.2 exact solution in the adiabatic limit 

As an approximation to H c , we truncate the chain at the 0-th 
site and obtain Ho(rj). It is noted that although fo tends to zero 
in the limit rj to 0, this truncation is not exact even in this limit. 
The reason is related to the singular coefficient in the operator 
X in Ho, to be discussed later. In order to study Ho(t]), we use 
the adiabatic approximation which is exact in the limit r\ — > 0. 
We check our results using the exact diagonalization method. 

To do the adiabatic approximation, we carry out another 
RWT H = VHoV- 1 where V = exp[KA(tf - b)](T z . Choosing 
a suitable parameter k, we obtain 



s A 1 , + 

Ho = _ 5 °~ x + 2 z + + 



(13) 



where A = Jr]o/n(l + so/rj). From Eq.dTTT), one obtains A ~ 
^(i-s)/2 an( j tne Sq ^ ^i-s j t j s eaS y t0 s jjow that A/so — > °° 
and A/ sq — > oo, meaning that the limit of rj = is the adiabatic 
limit. 

Using the replacement q - (/?' + b)/ V2 and p = i(b^ - 
we get 



with 



H = H a [q, <t z , o-J + H d [p], 



H a = -<x z - -cr x + — —Act z q + -s q 
2 2 V2 2 

1 2 
H d = -sop ■ 



(14) 



(15) 



In the adiabatic limit, the dynamical part H d does not play any 
role and the partition function of Ho and the average of a phys- 
ical observable O can be evaluated exactly by 



dqTre-P""^'-^, 

OO 



(O) 



i r°° 

mL^ ^^' (16) 



Here the trace is for the spin degrees of freedom 11271 . After 
redefining y/2Aq — > q, H a becomes 



H a — —cr z <t x H — o~ z q H — q 



(17) 



Here c = so/ A 2 . The rescaling of q changes a factor in Z and 
does not influence the physical results. In the limit r\ — > 0, one 
can prove using Eq.dTOt-dTTb that c = s/(2aa> c ) is a constant, 
independent of A. 



3 Results and Discussions 

3.1 magnetization <<x,) 

For physical analysis, we calculate the average of & z at zero 
temperature, the static susceptibility x(T) and the spin-spin cor- 
relation function C(w) at finite temperatures. Here we present 
the main results and put the details of derivation in Appendix 
A. 

The zero temperature magnetization for a > a c is deter- 
mined by the following equations, 

(<r z ) = -cq , (18) 

with qo being one of the solutions of the equation below, 

e + q 



cq = 



(19) 



V^ 2 + (e + q) 2 

From this expression, the critical coupling strength is obtained 
in the limit rj = as a c — As/(2o) c ), which is exactly the mean- 
field result ||23| . Here, due to the error introduced in truncating 
H c to Ho, Ho only qualitatively describe the effect of the bath. 
The critical behavior of (<r z ) can be identified exactly from 
Eq.(fT9b and it is also same as the mean-field result, 

(<r z Xe = 0)K(a-a c ) 1/2 (a > a c ); 

{a- z ){a = a c ) oc e l/ \ (20) 

This shows that our single-mode approximation gives the clas- 
sical critical exponent /3 = 1 /2 and 5-3, which are the correct 
result for the SBM in the regime < s < 1 /2. 

In Fig. 3, m = (<x,)/2 as functions of e are plotted on log 
scale for various a values at s — 0.3. For a < a c , m oc ^e. 
The susceptibility^ oc (a - a c ) 7 when a approaches a c . For 
a = a c , m oc e 1 ^ . For a > a c , m oc e°(a - a c *f in the small 
e limit. Numerical fit of data in Fig. 3 confirms the classical 
exponents f3 = 1/2, 6 = 3, and y — 1 (inset of Fig. 3). We 
observe that the m — e curves for a > a c and for a < a c always 
merge with the curve for a = a c in the regime e > e cr . Here e cr 
is the crossover value separating the quantum critical regime 
and the delocalized phase (for a < a c ) as well as the localized 
phase (for a > a c ). This feature implies the fulfillment of the 
scaling relation among/?, 5, and y. Indeed, using {a-a c y~> ' e cr = 
e® r (a - a c ~f = ej r 6 , we get the scaling relation 

r 



6-1 



(21) 



For the SBM, the same scaling relation has been confirmed for 
the regime < s < 1/2 [23|. The ED results are shown for 
several a < a c . Using r\ = 10~ 7 and Nb = 3500, we are able to 
reproduce the results of adiabatic limit at rj = only for small 
a values. For larger a, the error due to finite rj and Nb becomes 
large, and the exact results for rj = and Nb = oo are more and 
more difficult to obtain using ED. 
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Fig. 3. Order parameter m = (o"-)/2 at T - as functions of e for var- 
ious a's. The solid lines are results of adiabatic approximation. Solids 
dots are ED results using r\ = 1CT 7 and A^o = 3500. From top to bot- 
tom, a = 0.05,0.035,0.031,0.0302,0.03,0.0299,0.029,0.025,0.01. 
Among them a c = 0.03 has slope 1/3. Other parameters are s = 0.3, 
A = 0.2, and A = 2.0. Inset: m/e (for a < a c ) and m (for a > a c ) in 
the small e limit, as functions of \a - a c \. The dashed lines have slopes 
of 1/2 (for m) and -1 (for m/e). 



3.2 magnetic susceptibility v 

The magnetic susceptibility x(T) from the adiabatic approxi- 
mation reads 



- ( X 

'■£[HHf) + 7-(f)]'*'* 

B = 2 £^ sinh ^y-w dq . (22) 



Here u = ^A 2 + (e + q) 2 . The partition function Z is given by 
Z = 2 J" cosh(y)e-4«r^. (23) 

In Fig.4, ^(T) curves are plotted for different a values. We 
find similar structure as in Fig. 3. The ED results using 77 = 10" 7 
and Nb = 3500 agree with the adiabatic approximation in the 
low and high temperature regimes for small a. For larger ar's 
and in the crossover regime, obvious deviations occur due to 
finite 77 and Nb- In the low temperature limit, for a < a c , x x 
(a - aY y T°. Numerical fit confirms y — 1 (inset of Fig. 4). At 
a = a c , it is found that % oc T~ x with the exponent x = 1 /2. 
This is again same as the expected behavior for the SBM in 
the regime < s < 1/2. The BNRG calculation for the SBM 
gives x = s in this regime, due to the mass flow error. This error 
has been discussed and partly remedied by a modified BNRG 
method QjQ . For a > a c , X K ( a ~ a c)T~ l is consistent with 
X °c (o- z ) 2 /T in the localized phase and (3 = 1/2 (inset of Fig.4). 

We denote T* as the crossover temperature separating the 
quantum critical regime and the delocalized phase (for a < a c ) 
or the localized phase (for a > a c ). Similar observations of the 
scaling form in^(T) curves imply the following relations about 



Fig. 4. Magnetic susceptibility,^ = -3((<x z )/2)/<9e| E= o as functions 
of T for various ar's. The solid lines are results of adiabatic approx- 
imation. Solid dots are ED results. From top to bottom, a values are 
same as in Fig. 3. For a = 0.03, the slope is -1/2. Other parameters 
are 5 = 0.3, A = 0.2, e = 0, and A = 2.0. Inset: x (for a < a c) and T\ 
(for a > a c ) in the small T limit, as functions of \a — a c \. The dashed 
lines have slopes of -1 (for^) and 1 (for T%). 



T oc (a - a c y /x , 



(24) 



Since T* vanishes in the limit a — > a c as T* oc (or — ar f ) 1 '- v * [16], 
the above equations thus connect the critical exponent y, ft, x 
and y* by 

x 1 - x 



y t = - = 
r 



(25) 



Our results in Fig.4 and numerical fittings therefore confirm the 
classical exponents y* = 1 /2. 

It is noted that the results x(T) are different from that of 
the mean-field approximation, in which the Z2 symmetry of the 
SBM at e = is broken artificially. At a — a c , the mean-field 
approximation gives an exponential form for %\T), due to the 
finite spin gap induced by Acr x term. Here, although the spin 
gap is finite for a fixed adiabatic mode q, the integral over q 
effectively reduce the gap to zero at a = a c , leading to a power 
law in^(T). 



3.3 spin-spin correlation function C(co) 

For the spin-spin correlation function C(a>) defined below, 



C(t) = ^({cr z (t),cr z (0)}}H, 



1 r 

2n J_ a 



e l ""C(t)dt, 



(26) 



its expression in the adiabatic limit is obtained as 

C(cj) = w 6(cj) + C + (a>)0(cJ -A) + C_(<u)0(-w - A). (27) 
Here 

1 A 2 I _fc„2 _££„2\ . (PU>\ 



C + (a>) 



z V2|w| Vw 2 - A 2 



( e -^Ue-^)cosh(^). 



(28) 



6 F.-R. Liu, N.-H. Tong: Single Mode Approximation for sub-Ohmic Spin-Boson Model: Adiabatic Limit and Critical Properties 



(a) 

i 10'' 

- ~ ====:=:! ^~. 1 ° ! 


10" 10" 3 10" 2 ' 







io 8 10 5 io" 10 3 io z io' io° 



(B-A 





a=0.02 " 




a=0.035 ; 


■ (b) 




■ (c) 






10" 6 10" 5 


10 4 10" 3 10" ! 10" 


10 s 10" 5 


10" 10" 3 10" 2 


10"' 



<o-A (0-A 

Fig. 5. (a) C+(a>) versus u> - A for different a's at T = 10 5 . 
From top to bottom on the left: a = 0.01,0.025,0.03,0.031, and 
0.0315. The dashed lines are for a = 0.35,0.4, and 0.5 from left 
to right, respectively. Inset: wo (the weight of S(a>)) versus a - a c 
at T = 10~ 7 . The dashed line has a slope 1.0. (b) C + (co) versus 
ii) - A for a = 0.02 < a c at different T's. From top to bottom on 
the left: T = 10~ 6 , 10~ 5 , 10~ 4 , 10~ 3 , 10~ 2 . (c) C+(co) versus a) - A 
for a = 0.035 > a c at various T's. From top to bottom on the left: 
T = lQ-\3 x 10~ 4 ,2 x 10~ 4 , 10~ 4 . The dashed line is for T = 10~ 5 . 
Other parameters are 5 = 0.3, A = 0.2, e = 0, A = 2.0. 

and C-(oj) = C+(-u>). The function 6{co - A) is the step func- 
tion. q\ and are given by 

q\ = Vw 2 - A 2 - s, 
q 2 = - Vw 2 - A 2 - s. (29) 

It is seen that C(to) has a gap at \co\ < A except a delta peak at 
lo = 0. The weight of the delta peak Wo is 

W0 = l£ cosh (y)^ e " !e?2 ^ (30) 

Here u - ~\JA 2 + (e + q) 2 . From the definition, the sum rule 

XOO 
^ C(Lo)da) = 1 should hold. Using numerical integral, we 

have checked that the sum rule is fulfilled by Eq.(l27li-(l30il. 

The evolutions of C + (cj) with a or T are shown in Fig. 5. 

The weight Wo is finite at T > 0. At T — 0, it scales like (cr z ) 2 : 

wqoc (a - a c ) 213 for a > a c and wq — for a < a c , as shown in 

the inset of Fig. 5(a). At finite T, for a < a c , C+icS) divergence 

at the gap edge u> — A as {w - AY 1 ^ 2 . As a increases, this peak 

broadened until at a = a c , a Gaussian peak emerges. As a 

increases further, the weight continuously shifts from the gap 

edge to the Gaussian peak centered at to — (a/a c )A. These 

are shown in Fig. 5(a). For a fixed a, when T decreases, the 

weight of C + (cS) is concentrated either to the gap edge a> — A 

with square root divergence, or to the Gaussian peak at a> = 

(a/a c )A, depending on whether a < a c or a > a c , as shown in 

Fig. 5(b) and (c). Either peak becomes sharper with decreasing 



T. Eventually, at T — 0, it is expected that C+(cS) — 1 /26(a)- A) 
for a < a c , and C+{u>) = [(1 - wo)/2]5{a> - a/a c A) for a > a c . 

It is seen that C{a>) of the single-mode Hamitonian H a is 
nontrivial. At T = 0, it describes coherent quantum oscillations 
with frequency A in the delocalized phase, and with frequency 
{a/a c )A in the localized phase. At T > 0, the temperature de- 
pendent incoherent effects are reflected in the broadening of 
the delta peaks. These behaviors differ dramatically from that 
of the SBM. For the latter, at T = 0, C(u>) oc u> s in the small to 
limit. At the critical point, C(ai)(T = 0) oc <^r 5 . Using the RWT- 
based methods, dynamical properties of the SBM have been 
studied and rich phenomena disclosed 11 1II12L Recently, non- 
equilibrium dynamics was studied using the non-equilibrium 
BNRG QUI . The finite temperature dynamical properties of the 
sub-Ohmic spin-boson model is still a challenge to BNRG. 



3.4 discussions 

3.4.1 scaling relation 

Combining the information in Fig. 3 and Fig. 4, we conclude 
that these results are consistent with the scaling form of the 
critical part of the free energy F cri (j, e, T) ll3D (t = a - a c ), 

F crt (T,e,T) = Tf(— ^— ,— i-V (31) 

Here a T , aj, and a e are scaling constants for r, T, and e, respec- 
tively (32]]. Indeed, Fig. 6 demonstrates that Eq.PTTl is fulfilled. 
From Eq.OTT). all the critical exponents discussed above can be 
expressed using the scaling constants as 

/3 — (ot — a e ) /a T , 6 — a t j (aj — a e ) 

y — (2a e — aj) /a T , x — (2a e — aj) /aj, 

y* = a T /a T . (32) 

The scaling relations in Eq.(l2TTi and d25l ) can be regarded as 
consequences of Eq.PTTl. The critical properties of the single- 
mode Hamiltonian Ho(t] — > 0) can be summarized as a T /aj = 
1/2 and a e jaj - 3/4. 



3.4.2 problem in the BNRG crossover temperature T* 

Here we discuss the BNRG result of the the crossover tem- 
perature T* for the SBM. As implied in Eq.(l3ll T* oc T 1/r * 
|[T6ll is a temperature scale which marks the crossover between 
two distinct behaviors in F cr i{T). In BNRG, T* is obtained 
from the crossover point N cr in the flow of the energy lev- 
els using T* = A~ Ncr fl4l l 1411 151 . It has the universal meaning 
as in Eq.OTb. For the regime < s < 1/2, the validity of 
EqQTTl has been confirmed directly by QMC calculation Ifl6ll 
and indirectly by other studies, with the classical exponents 
a T jaj = 1/2 and 0^/07- = 3/4. This implies that T* oc r 2 
should hold inO < s < 1/2. However, previous NRG stud- 
ies gives T* cc t'/ s . We therefore conclude that besides the 
known defect of incorrect exponent /3, 6, and x in the regime 
< s < 1 12, a more direct defect in previous BNRG is the 
wrong crossover energy scale in the flow diagram. We have 
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Fig. 6. (a) \F(r,T,e = 0) - F(0,r,0)|/T versus \t\/T 1/2 for various 
T's in the localized and delocalized phases, (b) \F(t = 0, T, e) - 
F(0, T,0)\/T versus \e\/T 3/4 . Here t = a - a c . Other parameters are 
S =0.3, A = 0.2, A = 2.0 



checked that the boson state truncation does not influence the 
flow diagram [23 1. This problem in the crossover temperature 
T* is probably connected to the mass flow error IfTSll . Given 
that this issue is under active debate Il33ll . it is still a challenge 
for us to fully understand the problem with T* and remedy it 
within BNRG. 

The correct result v = 1 fs for < s < 1/2 obtained in pre- 
vious BNRG studies [4,. 14) is a result of cancelation of errors. 
It is obtained from the incorrect exponent y* = s fitted from the 
BNRG data, and the relation y* = 1/v which is not appropriate 
for systems above the upper critical dimension, v = 1 /s should 
be correctly obtained using y* = 1/2 and y* = 1/v + 1/2 — s, 
instead. 



3.4.3 truncation error in H (t] -> 0) 

The single-mode Hamiltonian Ho(rj — > 0) fails to quantitatively 
produce the correct a c for the spin-boson model. We would 
ask why is the truncation approximation invalid for H c , despite 
fo — > in the limit rj — > 0. To make this clear, we carry out the 
same unitary transformation V used for Eq.(fT3l) to the full chain 
Hamiltonian H c before it is truncated to Hq, H c = VH c V~ l . We 
obtain 



He = ^cr z - -cr x + - J — (1 + —)(b + bl)cr z 
2 V 7i rj u 



J][£nblb n + t n (blb n+i + bl +i b„)] 



n=0 



f o ho,, ,-k 



(33) 



In the limit 77 — > 0, the coefficient of the last term is nonzero, 
to/(2rf) y/rjo/n — > const. + 0. As a result, dropping out the 
term to(bTb\ +b' l bo) in H c means that the last term in the above 
equation is neglected, which is a strong approximation. 

H c is obtained from RWT and the subsequent inverse RWT 
It should be on the same boson basis as the original Hamilto- 
nian if the two RWT's cancel each other. Indeed, in the limit 
rj — > 0, all the terms containing bo or b Q disappears and <x- is 

coupled to (b\ + b x ) instead. The 0-th mode becomes a redun- 
dant degrees of freedom and H c recovers the standard NRG 
chain Hamiltonian. Numerical calculations confirms the fol- 
lowing relation between the coefficients, 



lim 


to 




= lim 


^0 






l]—>oo 


lim 


e« = 


c lim e n -\, 


17-0 






X} 


limf„ = 


c lim f„_i. 


,-»<> 




>!-> 


OS 



(34) 



c and c' tend to 1 in the limit A — > 1 . Therefore, in the limit 
A — > 1 where the discritization error disappears, H c return to 
the chain Hamiltonian H c exactly. 



3.4.4 possible NRG algorithm without boson state 
truncation error 

Since the implementation of BNRG based on H c does not pro- 
duce physical fixed point, here we propose a possible Wilson 
chain Hamiltonian which could avoid the boson state trunca- 
tion error, based on our present understanding. Starting from 
H c , we use the replacement q = (b^ + b)/ V2 and p = i(b^ — 
b)l V2. As done in Ref.[27 | for analysis of the adiabatic limit, 
we further carry out rescaling for p and q by defining 



P = 



V2 



—p 



1 = 




(35) 



In the limit 77 — > 0, according to Eq. dTTb . the coefficients of 
p 2 , {b\ + b'^p, and qcr z vanish. The remaining coefficients are 
nonzero constants. Then one gets a chain Hamiltonian without 
vanishing parameters, 

s A A 
H c (r) = 0) = -o- z - -cr x cosh (£) - -io- y sinh (£) 



e q 2 + t q (bi + b\j 

+00 

'Y J [e n b\b n + t n {b\b n+ \ + bl +l b n )]. 



(36) 



n=l 



Here,^ = — ip and the non-zero coefficients eo and to are 



to 



Urn 

r/^0 \ 7TT] Z 



i7-»0 \tj y n ) 



(37) 



The chain coefficients e„ = lim, ;= o e n {rj) and t n = lim 7= o t n (rj) 
are well defined quantities due to Eq.(l34ll. H c (rj = 0) has the 
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desired features: (i) in the limit A — 0, bosons are not dis- 
placed and the localized fixed point is correctly described with- 
out problem; and (ii) there is no vanishing parameters in the 
limit r] = as in H c . Therefore, it is a possible candidate for 
implementing BNRG without boson state truncation error. Fur- 
ther studies in this direction is in progress and the result will be 
discussed in a separate publication. 



4 Summary 

In this work, we explored the single-mode approximation for 
the SBM, using the RWT and the NRG-based transformations. 
Analytical results for the critical coupling strength a c , the mag- 
netic susceptibility x(T), and the spin-spin correlation function 
C(a>) at finite temperatures are obtained and confirmed by nu- 
merical results. We find that the derived single-mode Hamil- 
tonian gives mean-field a c and the critical exponents are clas- 
sical, /3 = 1/2, S = 3, y = \,x= 1/2, y* = 1/2. All these 
exponents are in agreement with those of the SBM in < s < 
1 /2. C(oj) has nontrivial behavior reflecting coherent oscilla- 
tion with temperature dependent damping effects due to envi- 
ronment. We discussed the problem in the crossover tempera- 
ture T* in the BNRG calculation, and propose a chain Hamilto- 
nian possibly suitable for implementing BNRG without boson 
state truncation error. 

This work is supported by National Program on Key Basic 
Research Project (973 Program) under grant 2012CB921704, 
and by the NSFC under grant number 1 1074302. 



A Appendix 

In this appendix, we present details for calculating (o%), x(T), 
and C(a>) using the adiabatic approximation which is exact in 
the adiabatic limit rj — > for Ho(t]). 

The matrix for H a {q) in Eq.dT7l under the o% eigen bases 
(IT), ID) reads, 



u - I 2 T 2 T 4 1 ' 2 
" ~ 1 A —£ — 1 -t- £-n 2 

' 2 2 2 4^ 



(38) 



The eigen energies are 



c u 

c u 
E 2 = -q 2 + -, 
C 2 



(39) 



with u = ^A 2 + (s + q) 2 . E\ = E g is the ground state energy. 
The corresponding eigen states are 



iPi = 



1 



V2m 



A [u + e + qY 1/2 
[u + £ + q] l/2 



and 



>l2u\-A[ 



[u + e + qf 2 
u + e + q] 



1/2 



(40) 



(41) 



The matrix elements of cr. are (¥'i|cr,|¥' 1 ) = -(e + q)/u, 
(¥ 2 \cr z \¥2) = -<5Pi|o-JPi>, {¥\<r z \<¥ 2 ) = <Fa|ff-«|Fi> = A/u. 
The order parameter (cr z ) is calculated using Eq.dToTi. At zero 
temperature, the expression reduces to 



<tr z > = <!P f 1 |tr z |!P r i>|^. 



9=90' 



(42) 



with go being the coordinate at the minimum of E g (q). qo can 
be solved by dE g {q)/dq = which leads to Eq.(fT9b. 

For susceptibility x(T) = -l/2<9(cr z )/<9e| e= o, one gets the 
expression 

dqTr(a z e-^) 

CO 



4Z 2 



e=0- 



(43) 



Using the matrix elements of cr z and E\ t 2, one gets the Eq.(l22l) 
in the main text. 

For the spin-spin correlation function C(a>), we start from 
the Lehmann expression 



C(u) 



i r°° 

— J dqj] [e~ PEm + e- pE ") Kn|^|m)| 2 5( W +£„-£ m ). 



(44) 



Using the cr, matrix elements and E\^, one gets 

(e + q?~ 



I r d q ( 



6(a)) 



(45) 



After the integral over q analytically, we get Eq.j27b-(f30b in the 
main text. 



B Appendix 

In this part, we present the numerical methods for solving the 
single-mode Hamiltonian Ho (Eq.lfTSl) in the adiabatic limit. In 
the limit rj tends to zero, A/eo ~ rf ( l ~ s » 2 — > oo. In this case, the 
boson occupation basis for Ho is not sufficient to produce an 
accurate result. We therefore carry out RWT to Ho and obtain 



6 t A y A -x 



(46) 



Here we use a — bo and a' = bl for convenience. cr ± = 



(<r x ± <T y )/2. x = a(a J - a) with a = A/e = y/T} /ji(l/eo + 
I/77). To diagonalize Ho, we need to first write down the ma- 
trix for the operator e ±x in the boson occupation basis {|f7j)} 
(m = 0, 1,2, ...,Nb - 1)- Then we construct Ho in the bases 
{|T>|m>,U>|m>Km = 0,l,...,Afc-l). 

Note that is the unitary operator that diagonalized the 
displaced harmonic oscillator h = a^a-a(a J +a), i.e.,e~ x he* = 
a* a— a 2 . Therefore we have (m|e*|n) = {m\n) a with \n) a the n-th 
eigen state of n, i.e., h\n) a — n- a 2 (m,n = 0,1, ...). 
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The eigen wave function of displace harmonic oscillator 
can be calculated using a number of methods. One could diag- 
onalize h in the occupation basis directly, or diagonalize^- first 
and then calculate its exponential. These methods are found 
to be inaccurate and slow. One could also use the recursive 
relation as done in Ref. [4|, or use the series summation as 
down in Ref. [6]. These methods are not satisfactory because 
for a > 10 1 and n > 10 2 , numerical unstability occurs. 

To overcome this problem, we resort to a new recursive al- 
gorithm for numerically calculating Laguerre polynomials. The 
eigen wave function of displaced harmonic oscilltor can be ex- 
pressed by Laguerre polynomial as [34] 



(m\n) a = (m\e a(a '- a) \n) (m,n = 0, 1,2, ...) 

'4z£ m -"V)a m -". 



or reduced by a fixed factor, say 1.02. Then the newly de- 
fined I^\x) is recalculated. This process iterates until \I ( „ z \x)\ 
is within the required range. Since P„ , (je) is small for a large x, 
this iteration is always able to find a suitable K„ (x) that makes 
within [l(T 200 , 10 2()0 ], avoiding the upper or lower over- 
flow. Finally, the matrix element of is calculated using the 
converged I^(x) and K% (x) as 



(m\n) a = /r«V) [P^V)] 1 " 4 "^ 



(52) 



Numerical test shows that this method can produce reliable re- 
sults for a ~ 10 3 and n ~ 10 5 without problem. 



( 4? ) References 



Here a — 1 is assumed even for a — 0. L„ (x) is the associ- 
ated Laguerre polynomial which obey the following recursive 
relations 



nl$\x) = (2n + z - 1 - ~{n + z~ 1)I£ 2 (jc), 

L ( f(x) =-x + z + l, 



L%\x) = L 



(48) 



This recursive relation is stable numerically, but L„ (a 2 ) over- 
flows for a > 50 and n > 500. At the same time, the prefactor 
e -i/2»- van i s hes for large a, leading to uncontrolled precision. 

To solve this problem, we design a self-adapting recursive 
algorithm to cancel the divergence. We define 



I%\x) = I%\x)[P®(x)] 
P z n (x) = e^¥' 2 



{n + z)\ 



(49) 



(z) 

The K„ (x) is a tunable parameter for each given x and n. The 
recursive relation for I ( n z> (x) reads 



,(Z) 



nij»(x) = (2» + z - 1 - jc)7«(jcH 1 



-(n + z- !)/»(*) 



n-2 v 



The initial conditions are 



[^iW] 



(50) 



/< z) (x) = (-x + z+l)[pf(x)] 



(51) 



To do the calculation, initially we set K^\x) = 1 /| In P [ f{x 



At the recursive calculation for some (n,z,x), if l/^WI > 10 



->200 

ic \n, 4., A,), 11 |J„ IV 

or (x)| < 10~ 200 occurs, A"„ (x) is respectively enlarged 
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